sf rastersf geometry typesst_point()st_linestring()st_polygon()st_multipoint()st_multilinestring()st_multipolygon()st_geometrycollection() # POINT (pt1 <- st_point(c(5, 2)))
# LINESTRING (ls1 <- st_linestring(rbind(c(5, 2), c(1, 3), c(3, 4))))
# POLYGON pl <- list(rbind(c(5, 2), c(1, 3), c(3, 4), c(5, 2))) (poly1 <- st_polygon(pl))
# MULTIPOINT mpmat <- rbind(c(5, 2), c(1, 3), c(3, 4)) (mp1 <- st_multipoint(mpmat))
# MULTILINESTRING
mll <- list(rbind(c(1, 5), c(4, 4), c(4, 1), c(2, 2)), rbind(c(1,
2), c(2, 4)))
(mls1 <- st_multilinestring((mll)))
# MULTIPOLYGON
mpyl <- list(list(rbind(c(1, 5), c(2, 2), c(4, 1), c(1, 5))),
list(rbind(c(0, 2), c(1, 2), c(1, 3), c(0, 3), c(0, 2))))
(mpoly1 <- st_multipolygon(mpyl))
# GEOMETRYCOLLECTION gcl <- list(ls1, pt1, mpoly1) (geomcol <- st_geometrycollection(gcl))
sfg to sfcsfgclass(pt1)
## [1] "XY" "POINT" "sfg"
class(mp1)
## [1] "XY" "MULTIPOINT" "sfg"
st_sfc()# sfc POINT pt1 <- st_point(c(5, 2)) pt2 <- st_point(c(1, 3)) (points_sfc <- st_sfc(pt1, pt2))
## Geometry set for 2 features ## geometry type: POINT ## dimension: XY ## bbox: xmin: 1 ymin: 2 xmax: 5 ymax: 3 ## epsg (SRID): NA ## proj4string: NA
class(points_sfc)
## [1] "sfc_POINT" "sfc"
is.list(points_sfc) # is this a list?
## [1] TRUE
is.list(mp1) # is this (multipoint) a list?
## [1] FALSE
st_geometry_type(points_sfc)
## [1] POINT POINT ## 18 Levels: GEOMETRY POINT LINESTRING POLYGON MULTIPOINT ... TRIANGLE
sfc to sfsfg and sfc only tell you about locations sf object is a data frame combines these two things:sfc) of this data frame lnd_point <- st_point(c(0.1, 51.5)) # sfg object
lnd_geom <- st_sfc(lnd_point, crs = 4326) # sfc object
lnd_attrib <- data.frame( # data.frame object
name = "London",
temperature = 25,
date = as.Date("2017-06-21")
)
lnd_sf <- st_sf(lnd_attrib, geometry = lnd_geom) # sf object
lnd_sf
## Simple feature collection with 1 feature and 3 fields ## geometry type: POINT ## dimension: XY ## bbox: xmin: 0.1 ymin: 51.5 xmax: 0.1 ymax: 51.5 ## epsg (SRID): 4326 ## proj4string: +proj=longlat +datum=WGS84 +no_defs ## name temperature date geometry ## 1 London 25 2017-06-21 POINT (0.1 51.5)
Let’s look at the CRS for the object we just created
st_crs(lnd_sf)
## Coordinate Reference System: ## EPSG: 4326 ## proj4string: "+proj=longlat +datum=WGS84 +no_defs"
CRS’s can be specified using an epsg code or a proj4string definition
epsg code refers to one well-defined coordinate reference system proj4string has various options and gives more flexibilityRead in a lat-long dataset
x_ll <- read.csv("data/confirmed-sightings.csv")
head(x_ll)
## Mountain Latitude Longitude Elevation ## 1 Jargalant Khairkhan 47.74961 92.48772 NA ## 2 Kharkhiraa-Turgen 49.76821 91.83882 1600 ## 3 Khasagt Khairkhan 46.77045 96.08086 2389 ## 4 Tsambagarav 48.72750 90.70322 2752 ## 5 Ikh Bogd NP 44.94847 100.03316 2224 ## 6 Altan Khukhii 48.69467 91.68228 2670
x_ll <- st_as_sf(x_ll, coords = c("Longitude", "Latitude"), crs = 4326)
st_crs(x_ll)
## Coordinate Reference System: ## EPSG: 4326 ## proj4string: "+proj=longlat +datum=WGS84 +no_defs"
d_ll <- st_distance(x_ll[1, ], x_ll[2, ])
x_utm <- st_transform(x_ll,
crs = "+proj=utm +zone=48 +datum=WGS84 +units=m +no_defs")
st_crs(x_utm)
## Coordinate Reference System: ## EPSG: 32648 ## proj4string: "+proj=utm +zone=48 +datum=WGS84 +units=m +no_defs"
d_utm <- st_distance(x_utm[1, ], x_utm[2, ]) c(d_ll, d_utm)
## Units: [m] ## [1] 229489.4 231902.2
occupancies <- st_read("data/Mongolia_SL.shp", quiet = TRUE)
head(occupancies)
## Simple feature collection with 6 features and 2 fields ## geometry type: POLYGON ## dimension: XY ## bbox: xmin: 26817.32 ymin: 5744980 xmax: 106817.3 ymax: 5804980 ## epsg (SRID): 32648 ## proj4string: +proj=utm +zone=48 +datum=WGS84 +units=m +no_defs ## Team Occ geometry ## 1 A 0.20980194 POLYGON ((66817.32 5784980,... ## 2 A 0.28545357 POLYGON ((86817.32 5784980,... ## 3 A 0.25328456 POLYGON ((66817.32 5764980,... ## 4 A 0.14932811 POLYGON ((86817.32 5764980,... ## 5 A 0.28437589 POLYGON ((26817.32 5744980,... ## 6 A 0.09210202 POLYGON ((46817.32 5744980,...
plot(st_geometry(occupancies))
plot(occupancies["Team"])
plot(occupancies["Occ"])
occupancies %>% st_crs() # check the CRS
## Coordinate Reference System: ## EPSG: 32648 ## proj4string: "+proj=utm +zone=48 +datum=WGS84 +units=m +no_defs"
occupancies <- st_read("data/Mongolia_SL.shp", quiet = TRUE) %>%
st_set_crs("+proj=utm +zone=48 +datum=WGS84 +units=m +no_defs")
occupancies %>% st_crs()
## Coordinate Reference System: ## EPSG: 32648 ## proj4string: "+proj=utm +zone=48 +datum=WGS84 +units=m +no_defs"
sights <- read.csv("data/confirmed-sightings.csv")
sights <- st_as_sf(sights, coords = c("Longitude", "Latitude"),
crs = 4326)
sights %>% st_crs()
## Coordinate Reference System: ## EPSG: 4326 ## proj4string: "+proj=longlat +datum=WGS84 +no_defs"
# reprojecting my_crs <- "+proj=utm +zone=48 +datum=WGS84 +units=m +no_defs" sights <- st_transform(sights, crs = my_crs)
plot(sights["Elevation"])
raster_filepath <- system.file("raster/srtm.tif", package = "spDataLarge")
new_raster <- raster(raster_filepath)
plot(new_raster)
mydata <- expand.grid(x = 1:100, y = 1:100) %>% mutate(z1 = x + y, z2 = x - y) myraster <- rasterFromXYZ(mydata) plot(myraster[["z2"]])
raster1 <- raster(nrows = 10, ncols = 10, res = 1, xmn = 0, xmx = 10,
ymn = 0, ymx = 10, vals = 1:100)
plot(raster1)
raster2 <- raster(nrows = 10, ncols = 10, res = 1, xmn = 0, xmx = 10,
ymn = 0, ymx = 10, vals = runif(100))
myrb <- brick(raster1, raster2)
plot(myrb)
dplyr verbs work with sf objects in intuitive ways
occ_gt90 <- occupancies %>% filter(Occ > 0.6) plot(occ_gt90)
Some specialised spatial functions like st_intersects
sights_per_cell <- st_intersects(x = occupancies, y = sights) sel_logical <- lengths(sights_per_cell) > 0 occ_with_sights <- occupancies[sel_logical, ] plot(occ_with_sights["Occ"])
occupancies %>% group_by(Team) %>% summarize(meanOcc = mean(Occ)) %>%
head()
## Simple feature collection with 6 features and 2 fields ## geometry type: GEOMETRY ## dimension: XY ## bbox: xmin: -773182.7 ymin: 4924980 xmax: 426827.5 ymax: 5804980 ## epsg (SRID): 32648 ## proj4string: +proj=utm +zone=48 +datum=WGS84 +units=m +no_defs ## # A tibble: 6 x 3 ## Team meanOcc geometry ## <fct> <dbl> <GEOMETRY [m]> ## 1 A 0.235 MULTIPOLYGON (((-53182.68 5484980, -53182.68 5504980, -33182.68… ## 2 B 0.473 MULTIPOLYGON (((-413182.7 5424980, -413182.7 5444980, -393182.7… ## 3 C 0.330 POLYGON ((-553182.7 5344980, -553182.7 5364980, -573182.7 53649… ## 4 D 0.407 MULTIPOLYGON (((-253182.7 5324980, -253182.7 5344980, -253182.7… ## 5 E 0.599 MULTIPOLYGON (((-473182.7 5104980, -473182.7 5124980, -453182.7… ## 6 F 0.318 MULTIPOLYGON (((386819.9 5144974, 386819.9 5164970, 366827.5 51…
Can also summarize geometries
occ_per_team <- occupancies %>% group_by(Team) %>%
summarize(meanOcc = mean(Occ),
geometry = st_union(geometry))
plot(occ_per_team[1:12,"meanOcc"])
sf, sp or raster.map0 <- tm_shape(occupancies) + tm_fill() map1 <- tm_shape(occupancies) + tm_borders() tmap_arrange(map0, map1, nrow = 1)
map2 <- tm_shape(occupancies) + tm_polygons()
map3 <- tm_shape(occupancies) + tm_polygons("Occ")
tmap_arrange(map2, map3, nrow = 1)
tm_shape(occupancies) + tm_polygons("Occ") + tm_text("Team",
size = 0.4)
tm_shape(occupancies) + tm_polygons("Occ") + tm_shape(sights) +
tm_dots("Elevation", size = 0.5)
Color settings are especially important for maps. Here use of the same palette obscures detail.
tm_shape(occupancies) + tm_polygons("Occ") + tm_shape(sights) +
tm_dots("Elevation", size = 0.5, palette = "BuGn")
tm_shape(occupancies) + tm_polygons("Occ") + tm_shape(sights) +
tm_dots("Mountain", size = 0.5, palette = "Set3")
tm_shape(occupancies) + tm_polygons("Occ") + tm_shape(sights) +
tm_dots("Elevation", size = 0.5, palette = "PRGn")
tm_shape(occupancies) +
tm_polygons(c("Occ", "Team")) +
tm_facets(ncol = 2)
geom_sf plots sf objectsggplot() + geom_sf(data = occupancies)
geom_sf plots sf objectsggplot() + geom_sf(data = occupancies, aes(fill = Occ))
sf objectsggplot() + geom_sf(data = occupancies, aes(fill = Occ)) + geom_sf(data = sights, colour = "red")
scale_XXX_brewer for discrete scales, scale_XXX_distiller for continuous scalesggplot() + geom_sf(data = occupancies, aes(fill = Occ)) + geom_sf(data = sights, aes(colour = Elevation)) + scale_colour_distiller(palette = "PRGn")
leafletleaflet() makes a map object that can be piped to other leaflet functionsleaflet() %>% addTiles()
leaflet() %>% addTiles() %>% setView(lng = 94.1, lat = 47.4,
zoom = 7)
leaflet() %>% addTiles() %>% setView(lng = 94.1, lat = 47.4, zoom = 4)
addProviderTiles()leaflet() %>% addProviderTiles("Esri.WorldStreetMap") %>%
setView(lng = 94.1, lat = 47.4, zoom = 4)
addProviderTiles()leaflet() %>% addProviderTiles("Esri.WorldImagery") %>%
setView(lng = 94.1, lat = 47.4, zoom = 4)
sf objects can be passed as additional map layersleaflet(data = st_transform(occupancies, crs = 4326)) %>%
addProviderTiles("Esri.WorldImagery") %>%
addPolygons()
leaflet(data = st_transform(occupancies, crs = 4326)) %>%
addProviderTiles("Esri.WorldImagery") %>%
addPolygons() %>%
addMarkers(data = st_transform(sights, crs = 4326))
pal <- colorNumeric("YlOrRd", domain = occupancies$Occ)
leaflet(st_transform(occupancies, crs = 4326)) %>%
addProviderTiles("Esri.WorldImagery") %>%
addMarkers(data = st_transform(sights, crs = 4326)) %>%
addPolygons(color = ~pal(Occ), fillOpacity = 0.4) %>%
addLegend(pal = pal, values = ~Occ) %>%
addMiniMap(position = "bottomleft")